Asymptotic behaviour of the Rayleigh— Taylor instabihty 
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We investigate long time numerical simulations of the inviscid Rayleigh- Taylor instability at At- 
wood number one using a boundary integral method. We are able to attain the asymptotic behavior 
for the spikes predicted by Clavin & Williams |15|| for which we give a simplified demonstration. In 
particular we observe that the spike's curvature evolves like while the overshoot in acceleration 
shows a good agreement with the suggested law. Moreover, we obtain consistent results for the 



prefactor coefficients of the asymptotic laws, 
interface profile near the spike. 



Eventually we exhibit the self-similar behavior of the 
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INTRODUCTION 

The Rayleigh- Taylor (RT) instability appears when, 
under gravity, an heavy liquid is placed over a lighter 
oneP|. This instability is crucial for our understanding 
of different phenomena in fluid mechanics: mixing, ther- 
mal convection and cited ref. herein) and also finger 
number selection in splashes Q. It is also important in 
inertial confinment fusion where the mass ablation pro- 
vides a stabilizing effect to the interface instability 0|. 
Without ablation, after the exponential growth of the 
perturbations due to the linear RT instability, nonlin- 
ear profiles develop through the formation of bubbles 
of lighter fluid rising into the heavier one and falling 
spikes of the heavier liquid penetrating the lighter one. 
In the general situations of viscous fluids which are im- 
miscible and/or have Atwood number not equal to unity 
{At = (ph- Pi)/ {Ph + Pi) with ph and pi being the den- 
sity of the heavier and lighter fluids respectively) , famous 
mushrooms- like structures grow for larger times 2, 5,Q. 
The limit of an inviscid fluid above a vacuum (At = 1) 
without surface tension plays a speciflc role since no sta- 
bilizing effects are present in the linear dynamics. Nu- 
merous theoretical and numerical work have focused on 
this idealized limit in order to track insights into the in- 
stability itself 7, 8, 9, 10, 11, 22J. It has been shown using 
a conformal mapping that a finite time singularity might 
appear in the conformal plane' 13^ and it is also suspected 
that for some sufficiently irregular initial conditions finite 
time singularities should also be observed in the physical 
plane. However, starting with sufficient ly s mooth initial 
conditions, the asymptotic dynamicsja llH u3 presents 
a constant velocity rising bubble separated by free falling 
tiny spikes as displayed on figure E] Although the rising 
bubble motion has been described using local properties 
of the flow[^, the asymptotic dynamics of the spikes 
is far from being well understood. The single mode ap- 
proach gives a fair description of the constant velocity of 
the rising bubble (uf, — \J gj (3/c) where g is the acceler- 
ation of the gravity and fc the wavenumber of the per- 



turbation) but gives only partial results for the spikejg. 
The fluid there obeys free fall dynamics to a good ap- 
proximation and the pressure fleld of the flow leads to an 
overshoot in the acceleration. The accelerated motion of 
the liquid stretches the spike geometry and one expects 
self-similar behaviour of the tip of the spikes. Recently, 




FIG. 1: Snapshots of the interface subject to the Rayleigh- 
Taylor instability for time ranging from f = to 10, start- 
ing with a small amplitude sine mode (left). On the right is 
shown the velocity of several points along the interface, non- 
dimensionalized with the stationnary bubble rising velocity 
^ g/ik, as a function of time. 

an asymptotic theory using a parallel flow description of 
the velocity field near the spikes has been constructed 
|l5|. The interface dynamics is nonlinear for large time 
and can be described using the theory of characteristics 
which gives rise to finite time singularity solutions. In 
the case of regular dynamics a self-similar description of 
the peak is obtained for large time: the maximal curva- 
ture of the interface at the peak tip is found to behave 
like the cubic power of time . Moreover, the spike po- 
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sition, following the free fall ^gt'^ at leading order, is 
shown to converge to the constant acceleration g with an 
overshoot in acceleration decreasing like t^^. In this let- 
ter, we present a numerical study of the Ray leigh- Taylor 
instability which focuses on the large time dynamics of 
the spikes in order to investigate the self similar dynam- 
ics predicted in 15]. We consider the dynamics for an 
inviscid liquid (heavy) with an exterior fluid of zero den- 
sity {At — 1) and no surface tension. The numerics use a 
boundary integral method (BIM later on) . Due to strong 
numerical instabilities, a careful treatment of the inter- 
face using conformal mapping is needed as explained be- 
low. The results are then shown and compared with the 
theory. 



ASYMPTOTIC ANALYSIS AND NUMERICAL 
METHOD 

We consider the two-dimensional motion of an inviscid 
fluid above a vacuum, subject to a negative acceleration 
—g. A periodic sine perturbation of the interface of wave 
number k is implemented as initial conditions. Neglect- 
ing surface tension, the equations of motion have no con- 
trol parameter after rescaling the time, the position and 
the velocity potential ip by factors \/gk, k and -y/ fc^/g 
respectively. The interface is described by y = a(x,t), 
where y is the direction along the gravity and x orthog- 
onal to it (see figure EJ. The velocity field U = {u,v) 
satisfies the dimensionless Euler equation 
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where P{x,y,t) is the pressure, Gy the non-dimensional 
acceleration due to gravity and the fluid density p = 1. 
The kinetic equation for the interface reads : 

da{x,t) da{x,t) 



dt 



dx 



with the velocity field (m, v) evaluated at the interface 
{x, a{x,t)). Starting at time t = with a small sine am- 
plitude interface, we observe for large time that the fluid 
particles located in the vicinity of the tiny spikes come 
from an almost free fall from the initial interface region. 
Therefore, following 0, we assume quasi-parallel steady 
flow for the velocity field which gives then in the tip re- 
gion |u| <C \v\ and : 
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with y ~ for large time. Writing a perturbation 
expansion of the velocity field in the tip region \x\ <^ y, 
we in fact consider: 



v=^/2{y + f{x,y,t)) 



with f{x,y,t) <C y. Taking a Taylor expansion in x of 
the perturbation /, we obtain by symmetry: 

,^^+Mgl + ^M0l + oix') 
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We limit our expansion to the second order in x for the 
velocity field later on. Incompressibility gives : 



1 , a(/o(y,o/v^) 

dy 



x + 0{x^). 



At the leading order (where we neglect even the pertur- 
bation f(x,y,t)) we obtain the following expression for 
the interface location: 



da{x,t) X da{x,t) 

dt ^2a{x,t) dx 



which can be solved using the methods of charasteristics 
(see ^3). Writing a{x,t) — t(| — j{x,t)) and noting 
that 7(0;, t) ^ t/2 in the spike region, we obtain, after 
linearisation : 



d^{x,t) xd^{x,t) 



dt 



t dx 



= 



which has self-similar solution of the form ^{x. t) — 9{xt). 
A first conclusion can be drawn about the curvature of 
the interface at the tip, k = —d'^a/dx'^\x=vn which is thus 
found to increase as the cubic power of time : 



(1) 



The next order terms of the expansion allow the de- 
termination of the function foiy,t) near the tip. Using 
the constant value of the pressure at the interface we use 
the projection of the Euler equation at the interface on 
its local tangent : 

du ^ da{x,t) dv da{x,t) 
dt dx dt dx 

Since on the interface dP{x,a{x,t),t)/dx = 0. We 
develop this equation at first non-zero order (which 
will end up to be the first order in x) with the ex- 
pansion e{xt) = 61(0) -I- xH^9" {0)/2 + 0{x*). Re- 
membering that I/I <C y, we can neglect also 
the larsge scale terms d'^{fo{y,t)/y/2y)/dtdy and 
^/2yd'^{fo{y, t) / / dy"^ with respect to the others. We 
obtain finally for the tip position y = ys : 



dfo{ys,t) 

dt 



-dfo{Vs,t) _ dfo{ys,t) 



dy 



dt 
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Recalling that: 



tip acceleration at leading order: 

d'^Vs , 1 dfo{ys,t) 



■^s + ^^^^tII^ we obtain for the 



df^ 



dt 



1 + 



t^e"{o) 



(2) 
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which corresponds to an overshoot in the spike accelera- 
tion decreasing as the fifth power of time. 

The numerical method is elaborated using the incom- 
pressible and potential properties of the flow. The veloc- 
ity field can thus be evaluated everywhere when the veloc- 
ity potential is known on the interface thanks to Cauchy's 
theorem, in the spirit of pionnering works |l6lll7Lll8lll9l| . 
The non-dimensional Bernoulli equation on the free sur- 
face reads : 



dip 
'dt 



1 



(3) 



where the velocity potential 95 is a harmonic function in 
the fluid domain 17 : 



Av? = 



(4) 



The kinematic condition on the free surface expresses 
the fact that fluid particles move with the same normal 
velocity than the free surface itself : 



Vip ■ n 



(5) 



Knowing ip on the free surface at a given time-step, we 
search for the solution of equation Q that satisfies this 
boundary condition (|3J). We use the complex potential 
/3(z) — ip + lip and the conformal map f{z) — exp{—iz) 
(Cf. Figure I2J), where z = x + iy and ip is the stream 
function. The conformal map transforms the periodic 
domain il into the closed domain AI. Since ip is harmonic 
inside fl, f3{z) is analytic inside fl and therefore 7(C) = 
f3{f{z)) is analytic inside M. Using Cauchy's theorem, 
we obtain a Fredholm equation of the second kind for the 
stream function ijj which is solved using discretization of 
the free surface {dil and thus dM). This linear system 
of equations is solved using a LU decomposition. Once 
we know ip on each point on dM, the complex velocity 
of each marker in the physical plane is given by : 



dp 
dz 



(6) 



where u and v are the horizontal and vertical velocities 
respectively. This complex velocity is computed with a 
finite difference scheme using the values of the complex 
potential on the collocation points on 911. The posi- 
tion of the surface markers (kinematic condition) and the 
value of the velocity potential on each of these markers 
(Bernoulli equation) are then updated in time using a 
fourth order Runge-Kutta method. 




FIG. 2; Conformal map used to transform the physical peri- 
odic plane Q into a closed domain M. 



sine-mode. The unavoidable numerical noise cannot be 
damped by the numerics and the calculations always end 
up subject to numerical instabilities. Nevertheless, we 
emphasize that the numerical scheme used here is re- 
markably robust and can be accurately evolved to reach 
the large time where the scalings predicted by the theory 
|l5| are valid. Comparing our simulations with recent 
numerical works0 0, 13> ''^^ have been able to run the 
dynamics at least twice as far which corresponds roughly 
to an increase of a factor of 8 in the tip's curvature. 

The position of the spike is shown on figure |3| as 
function of time. We observe that the asymptotics 
dynamics are very well approximated by the relation 
Us = \g{t — io)^ as shown in the inset to the figure with 
to — 3.74. This remarkable behavior, in good agreement 
with the free fall hypothesis, suggests that Iq is the time 
delay accounting for the initial exponential development 
of the instability. We will therefore present further data 
on the curvature dependance and the acceleration of the 
tip as functions of this delayed time t ~ to instead of t. 
The curvature Kg at the tip is then shown on figure 01 




RESULTS AND DISCUSSIONS 

We have performed numerical simulations of the 
Rayleigh- Taylor instability using the numerical method 
described above. We start with a small amplitude 



FIG. 3: Position of the spike j/s(t) as a function of time. The 
inset shows in a log-log plot of the spike position (black curve) 
as function of time t — to with to = 3.74 obtained by a second 
order polynomial fit of j/s . the dashed line shows the expected 
behavior ^t'^. 
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The large time asymptotic behavior is similarly found to 
follow the cubic law (see equation with 6"{0) = 1.5. 
In addition, the acceleration of the tip is computed by 
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FIG. 4: Spike curvature Ka calculated at the tip y = j/s as 
function of the delayed time t — to in a log- log plot. The 
dashed line displays the cubic law Q with ^"(0) — 1.5. 

finite differences on the tip velocity and the overshoot in 
the acceleration is presented on figure|Sl We observe that 
the results look noisier than the two previous ones. Two 
factors can explain such noise: firstly, we are looking to 
a finite difference which decreases to zero so that the nu- 
merical errors are relatively more important. However, 
we note that the overshoot in acceleration shows a good 
agreement with the law, noting that no adjustable 
parameter is used in this comparison. Moreover, the self 

lo" [ ^ ^ ^ ^ ^ — ^ — ^ — ^ 




FIG. 5: Overshoot in acceleration, defined as the difference 
between the tip acceleration and the gravity. The plot is in 
log- log scale and with the delayed time t — to- The dashed line 
shows the theoretical prediction Q using the value of 9"{0) 
obtained from figure 2] 

similar structure of the interface near the tip has been 
exhibited on figureEl We observe after the proper rescal- 
ing on the left part of the figure that the interface profiles 
collapse onto a single curve near the spike. 




FIG. 6: Self-similar structure of the tip: the interface profile 
around the spike have been superimposed on the right side 
of the figure for different time t ranging from 4 to 12. The 
left side of the figure shows the same curves rescaled by factor 
l/(t — to) and (t — to) for the x and y coordinates respectively, 
following the scaling behavior predicted by the theory. 

We have thus exhibited large times numerical simu- 
lations of the Rayleigh- Taylor instability which present 
asymptotic scaling behavior in agreement with theoret- 
ical predictions using Taylor expansions of the free fall 
velocity field at the spike Although our numerics al- 
ways stops due to numerical instability, we have been able 
to reach large time enough to exhibit the cubic power in 
time dependance for the spike curvature and the inverse 
of the quintinc power of time decreasing of the overshoot 
in acceleration. 
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